Bose-Einstein Quantum Statistics and the Ground State of Solid 4 He 
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The ground state of solid 4 He is studied using the diffusion Monte Carlo method and a new 
trial wave function able to describe the supersolid. The new wave function is symmetric under the 
exchange of particles and reproduces the experimental equation of state. Results for the one-body 
density matrix show the existence of off-diagonal long-range order with a very small condensate 
fraction ~ 10 -4 . The superfluid density of the commensurate system is below our resolution thresh- 
old, p B /p < 10 -5 . With a 1% concentration of vacancies the superfluid density is manifestly larger, 
Ps/p = 3.2(1) ■ 10- 3 . 

PACS numbers: 67.80.-s,02.70.Ss,67.40.-w 
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Solid 4 He is far from being a classical crystal as proved 
by its high degree of anharmonicity, large kinetic energy 
per particle (T ~ 24 K), and significant displacements, 
measured by the Lindemann ratio 7 ~ 0.26 (to be com- 
pared with classical solids, with 7 ~ 0.03 near melting). 
The counterintuitive possibility of simultaneous solid or- 
der and superfluidity in solid 4 He has attracted the at- 
tention of both theory and experiment from long time, 
starting with the theoretical proposal of Andreev and Lif- 
shitz [l| about a possible supersolid phase related to the 
presence of a finite concentration of vacancies. Finite val- 
ues of the superfluid density and/or condensate fraction 
in solid 4 He would be a macroscopic effect induced by 
Bose-Einstein statistics. After several unfruitful experi- 
ments to detect superfluid signals in solid 4 He, Kim and 
Chan reported in 2004 the first evidence of non-classical 
rotational inertia (NCRI) both in a confined environ- 
ment [2] and in bulk From then on, several other 
experimental groups have observed NCRI using different 
samples containing small or ultra small 3 He concentra- 
tions, in a simple crystal or in a polycrystal, and using 
several annealing schemes There is an overall agree- 
ment of all the data concerning the onset temperature 
T = 75 — 150 mK at which the superfluid fraction be- 
comes zero, the lowest value corresponding to ultra pure 
samples (only 1 ppb 3 He). Much more controversial is 
the value of the superfluid density since the experimen- 
tal values reported so far change by more than one order 
of magnitude (fl s /p — 0.03 — 0.5%) depending on the pu- 
rity, annealing, conditions in which the crystal is grown, 
etc. This high dispersion has led to think that the super- 
fluid signal observed in solid 4 He is probably due to the 
presence of some defects in the crystal, which could be of 
different nature: dislocations, vacancies or grain bound- 
aries. In fact, superfluid flow has been detected when 
grain boundaries are present [5j and also recent measure- 
ments of the specific heat has shown a broad peak at the 
same onset temperature Tq Q. 



Theoretical calculations at very low temperatures, 
based on the path integral Monte Carlo method (PIMC), 
have not been able to reproduce the experimental find- 
ings on the supersolid. PIMC results show that a com- 
mensurate perfect crystal does not exhibit neither su- 
perfluid fraction 0] nor condensate fraction Finite 
values of p s /p have been observed only when disorder is 
introduced in the form of a glassy phase [9( or through 
dislocation lines [lOj, but puzzlingly no signal of the su- 
perfluid transition is obtained over the temperature range 
observed in experiments. As the critical temperature 
of a supersolid is not known, conceptually any finite- 
temperature method can not provide the definite answer 
and one has to resort to a strictly zero-temperature cal- 
culation. In this letter, we propose a new trial wave func- 
tion which allows simultaneously for spatial solid order 
and Bose symmetry, and with the benefit of a simple use 
for importance sampling in diffusion Monte Carlo (DMC) 
calculations. 

Structure and energetic properties of solid 4 He at zero 
temperature have been widely studied in the past using 
the Nosanow-Jastrow trial wave function, 
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N being the number of particles and lattice sites (com- 
mensurate crystal), f(r) a two-body correlation function, 
and g(r) a one-body localization factor which links every 
particle i to its site I. The wave function \1/nj leads to an 
excellent description of the equation of state and struc- 
ture properties of the solid [llj, but it can not be used 
for estimating properties which depend directly on the 
Bose-Einstein statistics since it is not symmetric under 
the exchange of particles. The symmetry requirement 
can be formally written as 
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with a sum over all the permutations P(J) of lattice 
sites. This wave function has been used in the past in 
variational Monte Carlo (VMC) calculations by introduc- 
ing an explicit sampling over the permutation space [l2T ] . 
However, this sampling is inherently inefficient and diffi- 
cult to incorporate in a DMC code. To avoid the com- 
plexity of random walks in the permutation space one 
can approximate ^pnj by flii ] 
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This wave function fulfills the symmetry requirement but 
it is not well-suited for genuinely describing solid 4 He. 
The problem with using ^lnj as importance sampling 
function in DMC is that the solid phase melts, arriving 
to liquid (L) or glassy configurations. This outcome can 
be understood by analyzing the behavior of the localiza- 
tion factor (second term in (J3])) under the per-site occu- 
pancy number. Contrary to what is physically reason- 
able, multiple occupancy per site is only suppressed via 
the Jastrow part (first term in ([3])). To preserve the solid 
structure it is necessary to take into account, already in 
the localization factor, a penalization accounting for the 
voids in the original sites created by multiple occupancy, 
a feature missing in the ^lnj wave function. 

We introduce a new type of wave function vE'sn.i in 
which this feature is present and at the same time the 
necessary requirements of solid order and Bose symmetry 
are fulfilled. The key point is to use the site occupancy as 
the building block for the localization factor, thus voids 
are unequivocally taken into account. Our model takes 
the form 
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where the product in the second term runs over sites in- 
stead of particles, and the localization factor properly 
suppresses the voids arising due to double occupancy, 
a right behavior which is also present in the more gen- 
eral wave function ^pnj- There have been other pro- 
posals for the study of the supersolid which do not rely 
on the symmetrization of the Nosanow-Jastrow model: 
a Bloch-like function inspired in the band theory 
for electrons, and the shadow wave function The 
first model was used in the past in VMC simulations of 
Yukawa solids [13] and proved to be difficult to optimize 
and provided worse energetic results than the NJ model. 
More successful has been the shadow wave function which 
has been applied recently in VMC calculations of super- 
solid 4 He [14j . However, it has been never applied as im- 
portance sampling wave function in a DMC calculation. 
Recently, it has been implemented in the path integral 
ground state (PIGS) method to eliminate the variational 
constraints and applied to the study of two-dimensional 
solid 4 He [H. 
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FIG. 1: DMC results for the energy per particle of solid 4 He 
as a function of the density. Squares and diamonds stand for 
the symmetric (SNJ) and non-symmetric (NJ) wave functions, 
respectively. The solid line is a polynomial fit to the SNJ 
results. Experimental data from Ref. 17] is plotted with 
solid circles. 



The ground state of solid 4 He has been studied by 
applying the DMC method to N atoms in a simu- 
lation box with periodic boundary conditions. The 
DMC method solves stochastically the imaginary-time 
(r) Schrodinger equation, providing exact results for bo- 
son systems within controllable statistical errors. When 
t — > oo, the ground state dominates and one has a collec- 
tion of walkers Rj = {r 1; . . . , Tn} which follow the prob- 
ability distribution function (\E , 4'), ty and \P being the 
ground-state wave function and the trial wave function 
for importance sampling, respectively. The short-time 
Green's function according which the walkers evolve is 
accurate to order (At) 3 [l6| and internal parameters of 
the calculation such as the mean population of walkers 
and time step At have been adjusted to eliminate any 
bias. 

We have carried out DMC simulations of solid 4 He us- 
ing an hep lattice with N = 180 and 448 atoms and a fee 
one with TV = 108; the range of densities analyzed starts 
slightly below melting, p = 0.470 ct" 3 (a = 2.556 A) 
with pressure P — 21 bar, and ends at p = 0.600 cr~ 3 
with P = 160 bar. The Jastrow factor in <3>nj CQ) and 
^snj © is of McMillan type, f(r) = exp[-0.5(6/r) 5 ], 
and the Nosanow term is in both cases a Gaussian, 
g(r) — exp(— 0.5/3r 2 ). The parameters b and (3 have 
been optimized using VMC, the optimal values being 
b = 1.12 a and (3 = 7.5 <r~ 2 , and we have neglected 
their slight density dependence. 

In Fig. 1, we show DMC results for the energy per 
particle E/N as a function of the density for both the 
symmetric SNJ © and non-symmetric NJ ([I]) wave func- 
tions. The calculations used an hep lattice with iV = 180 
atoms and size corrections to the energy were estimated 
using the iV-dependence observed in VMC calculations of 
the same system. In a recent work llj, we have verified 



that this procedure improves the equation of state with 
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FIG. 2: Static structure function S'(fc) at density 
p = 0.491 a -3 (AT = 180) calculated with the right sym- 
metry (SNJ) (squares) and with the NJ model (triangles). 
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FIG. 3: One-body density matrix p(r)/p of solid 4 He at densi- 
ties p — 0.491 (circles) and p — 0.535 a" 3 (squares). Full tri- 
angles stand for a calculation at p = 0.491 <r -3 with N = 180 
particles and open triangles for the NJ result. 



respect to the standard approach of assuming that the 
system is homogeneous beyond L/2, with L the length of 
the simulation box. As is known, the NJ model is able to 
reproduce the experimental equation of state with high 
accuracy; our present results confirm this feature as can 
be observed in Fig. 1. More importantly, the results 
obtained with the right symmetrization (SNJ) are statis- 
tically indistinguishable of the NJ ones and arc in agree- 
ment with experiment. On the contrary, we have carried 
out some simulations with the LNJ model ([3]) and found 
that the energies are larger than the SNJ and NJ ones 
and the solid order is lost, the spatial structure resem- 
bling a glassy state. 

One of the most clear signatures of the solid phase 
is drawn from the form of the static structure function 
S(k) = l/iV(p k p_ k ), with p k = Etie 41 " 1 - In Fig. 2, 
we show DMC results of S(k) for hep solid 4 He calculated 
using the pure estimator method. The results presented 
correspond to a density p = 0.491 a~ 3 and have been ob- 
tained for the symmetric (SNJ) and non-symmetric (NJ) 
wave functions. Both results look rather similar, with 
the three main singular peaks corresponding to the re- 
ciprocal lattice points clearly visible. The height of these 
peaks for the symmetric system is slightly smaller than 
for the NJ one but these results confirm that the solid 
order is preserved by vE'sn.t- The small decrease in the 
strength of the main S(k) peaks is probably related to 
the atomic diffusion we have observed along the simu- 
lations. By monitoring the distance of any particle to 
its initial position along the simulation one can know if 
there is atomic diffusion in the crystal as a consequence 
of particle exchanges. Our results show that ~ 15% of 
the atoms present displacements r/a > 2, with a the lat- 
tice constant; obviously, the same measurement for the 
NJ case does not show any of them. On the other hand, 
the relevance of these exchanges is not observed in the 
density profile around the lattice points /i(r). The results 
of u(r) obtained with SNJ and NJ are statistically indis- 



tinguishable and therefore both Lindcmann ratios have a 
common value 7 = 0.26, in agreement with experimental 
data. 

A well-known drawback of the NJ model is the im- 
possibility of answering the fundamental question about 
the possible existence of off-diagonal long range order 
(ODLRO) in solid 4 He. Instead, the symmetrized model 
SNJ fulfills the right Bose-Einstein statistics and there- 
fore is able to provide this information. Quantitatively, 
ODLRO is measured by the condensate fraction no, that 
is estimated through the asymptotic behavior of the one- 
body density matrix p(r)/p, uq = lim r ^oo p(r)/p. Re- 
sults for p(r)/p at densities p = 0.491 and p = 0.535 a~ 3 , 
corresponding to pressures P = 31 and 68 bar respec- 
tively, are reported in Fig. 3. At both densities one 
can see unambiguously the existence of ODLRO, and 
from the asymptotic behavior we estimated the values 
of n : 4.3(2) • 10~ 4 and 1.7(4) • 10~ 4 at p = 0.491 and 
p = 0.535 c~ 3 , respectively (figures within parentheses 
correspond to the errors). For comparison, we show in 
the figure results for p(r) obtained with the NJ model. In 
order to analyze the influence of the number of particles 
used in the simulation, we have carried out calculations 
of p(r) with N = 180 and 448. The results obtained 
are plotted in Fig. 3 at the lower density; within the 
statistical error both estimations lead to the same value 
of no- All the results plotted in Fig. 3 have been ob- 
tained using extrapolated estimators fulfilling the con- 
dition that extrapolations of the same accuracy, namely 
n ~ 2nQ nix - n™ and n ~ (n nix ) 2 /n ' ar (mix and var 
stand for mixed and variational estimators, respectively), 
coincide. This equality is achieved by slight variation of 
the parameters of the trial wave function ((4]) which oth- 
erwise does not modify the energy results. It is worth 
mentioning that we have made some VMC calculations 
removing the Jastrow part of the SNJ trial wave function 
and, in all cases, a finite value of the condensate fraction 
is observed. This result points to a property which seems 
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FIG. 4: Diffusion of the center of mass aD a {r) as a function 
of the imaginary time r for solid 4 He at p — 0.491 cr~ 3 . The 
slope of the lines is directly the superfluid density p a j p. The 
solid line stands for the commensurate solid; the dotted line 
for the case with N — 1 sites; the dashed line for a solid with 
a ~ 1% concentration of vacancies. 

to be inherent to the symmetrization of the NJ model. 

From the different experimental measures carried out 
on the supersolid, it seems clear the the presence of va- 
cancies, dislocations or other defects can be crucial to 
understand the dispersion on the superfluidity results. 
We have used the SNJ model to determine this influence 
in the condensate fraction and superfluid density for two 
particular cases: a vacancy, i.e., a system with N lattice 
sites and N — 1 particles, and the absence of a lattice 
site, i.e., N particles and TV — 1 sites. Our simulations 
show that with ~ 1% vacancy concentration the conden- 
sate fraction raises nearly a factor two with respect to the 
commensurate case: at p — 0.491 cr~ 3 , n = 8.7(4) • 10~ 4 . 
On the contrary, the effect of removing one of the sites in 
a 1% fraction does not modify the commensurate value 
n = 4.3(2) • 1CT 4 . 

The superfluid density of a bosonic system can 
be calculated with DMC by extending the winding- 
number technique, originally developed for PIMC cal- 
culations, to zero temperature (HI- Explicitly, p s /p = 
linv-nxj q:(D s (t)/t), where a = N/6Dq with Dq = 
h 2 /2m, and D s (t) = ((R CM (t) - R C m(0)) 2 > with R CM 
the center of mass of the particles in the simulation box. 
In Fig. 4, results for the function aD s (r) at a density 
p = 0.491 (7~ 3 are plotted as a function of the imaginary 
time r. The slope of this function is directly the super- 
fluid density p s /p. Our results show that p s /p for the 
commensurate solid is smaller than our precision limit, 
that we have established at 1 • 10 -5 . We have carried 
out simulations with different number of particles with 
the same result. Also we have calculated D s (t) at higher 
pressure and at a pressure below melting where the solid 
is metastable obtaining identical conclusion. No signifi- 
cant difference with the regular system has been detected 
for the case with N — 1 sites, also shown in Fig. 4. On 
the other hand, when vacancies are present in the crys- 



tal with a ~ 1% concentration a superfluid signal clearly 
emerges; from the slope, p s /p = 3.2(1) ■ 10~ 3 . 

To summarize, we have introduced in the microscopic 
description of solid 4 He a new wave function (SNJ) which 
properly symmetrizes the well-known Nosanow-Jastrow 
model. Our results, based on the essentially exact DMC 
method, prove that the SNJ trial wave function used 
for importance sampling in the method is able to simul- 
taneously reproduce the experimental equation of state 
(as the NJ does) and predict results for the conden- 
sate fraction and superfluidity due to its Bose symmetry. 
The one-body density matrix shows ODLRO with a very 
small condensate fraction, uq — 4.3(2) • 10~ 4 near melt- 
ing. A finite but smaller value has been obtained by Galli 
et al. [LJ] using VMC with the shadow wave function. On 
the contrary, PIMC simulations at finite temperature do 
not observe ODLRO [8J]. An important conclusion of 
our work is that the zero-temperature upper-bound for 
the superfluid fraction p s /p < 1 ■ 10~ 5 19] is incom- 
patible with recent experimental measurements(3 • 10~ 4 - 
5-ht 3 ) i, ruling out an explanation of the superfluid 
signal as that of a supersolid in a perfect crystal. The 
introduction of 1% vacancies in the system increases the 
condensate fraction and a clear signature of superfluidity 
is detected. Work is in progress to understand micro- 
scopically the connection of other defects and/or disorder 
with both p s and no. 
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